Code
knitr::opts_chunk$set(cache = TRUE)Exploration of dataset
knitr::opts_chunk$set(cache = TRUE)Note: data last modified 7-Nov-2024 This dataset contains the scored moult data for all usable photos from the 2021-2022, 2022-2023, and 2023-2024 breeding seasons at Kaikoura
# 2021-2022 season data
ggplot(data = dat %>% filter(season == 2021 & sex == "F") %>% mutate(score = as.numeric(score),
rings_comb_ = reorder(rings_comb, n_photos, decreasing = TRUE))) +
geom_point(aes(y = 1, x = date, fill = score),
pch = 21, color = "black", size = 3) +
facet_wrap(. ~ rings_comb_, ncol = 1, strip.position = "right") +
scale_fill_gradient(high = "#cc4c02", low = "white") +
theme_bw() +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
strip.text.y.right = element_text(angle = 0)) +
scale_x_date(date_labels = "%W", expand = c(0.01, 0.01),
date_breaks = "3 week",
limits = c(as.Date("2021-07-05"), as.Date("2022-05-01"))) +
xlab("week") +
ggtitle("Females from the 2021-2022 breeding season in Kaikoura")ggplot(data = dat %>% filter(season == 2021 & sex == "M") %>% mutate(score = as.numeric(score))) +
geom_point(aes(y = 1, x = date, fill = score),
pch = 21, color = "black", size = 3) +
facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
scale_fill_gradient(high = "#cc4c02", low = "white") +
theme_bw() +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
strip.text.y.right = element_text(angle = 0)) +
scale_x_date(date_labels = "%W", expand = c(0.01, 0.01),
date_breaks = "3 week",
limits = c(as.Date("2021-07-05"), as.Date("2022-05-01"))) +
xlab("week") +
ggtitle("Males from the 2021-2022 breeding season in Kaikoura")# 2022-2023 season data (not as good coverage as the previous season)
ggplot(data = dat %>% filter(season == 2022 & sex == "F") %>% mutate(score = as.numeric(score))) +
geom_point(aes(y = 1, x = date, fill = score),
pch = 21, color = "black", size = 3) +
facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
scale_fill_gradient(high = "#cc4c02", low = "white") +
theme_bw() +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
strip.text.y.right = element_text(angle = 0)) +
scale_x_date(date_labels = "%W", expand = c(0.01, 0.01),
date_breaks = "3 week",
limits = c(as.Date("2022-07-05"), as.Date("2023-05-01"))) +
xlab("week") +
ggtitle("Females from the 2022-2023 breeding season in Kaikoura")ggplot(data = dat %>% filter(season == 2022 & sex == "M") %>% mutate(score = as.numeric(score))) +
geom_point(aes(y = 1, x = date, fill = score),
pch = 21, color = "black", size = 3) +
facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
scale_fill_gradient(high = "#cc4c02", low = "white") +
theme_bw() +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
strip.text.y.right = element_text(angle = 0)) +
scale_x_date(date_labels = "%W", expand = c(0.01, 0.01),
date_breaks = "3 week",
limits = c(as.Date("2022-07-05"), as.Date("2023-05-01"))) +
xlab("week") +
ggtitle("Males from the 2022-2023 breeding season in Kaikoura")# 2023-2024 season data
ggplot(data = dat %>% filter(season == 2023 & sex == "F") %>% mutate(score = as.numeric(score))) +
geom_point(aes(y = 1, x = date, fill = score),
pch = 21, color = "black", size = 3) +
facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
scale_fill_gradient(high = "#cc4c02", low = "white") +
theme_bw() +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
strip.text.y.right = element_text(angle = 0)) +
scale_x_date(date_labels = "%W", expand = c(0.01, 0.01),
date_breaks = "3 week",
limits = c(as.Date("2023-07-01"), as.Date("2024-05-01"))) +
xlab("week") +
ggtitle("Females from the 2023-2024 breeding season in Kaikoura")ggplot(data = dat %>% filter(season == 2023 & sex == "M") %>% mutate(score = as.numeric(score))) +
geom_point(aes(y = 1, x = date, fill = score),
pch = 21, color = "black", size = 3) +
facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
scale_fill_gradient(high = "#cc4c02", low = "white") +
theme_bw() +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
strip.text.y.right = element_text(angle = 0)) +
scale_x_date(date_labels = "%W", expand = c(0.01, 0.01),
date_breaks = "3 week",
limits = c(as.Date("2023-07-05"), as.Date("2024-05-01"))) +
xlab("week") +
ggtitle("Males from the 2023-2024 breeding season in Kaikoura")ind_breeding_scores %>%
group_by(rings_comb, sex) %>%
summarise(max_breeding_score = max(max_breeding_score)) %>%
ungroup() %>%
group_by(sex, max_breeding_score) %>%
summarise(n = n()) %>%
bind_rows(., data.frame(sex = c("F", "M"), max_breeding_score = c(7, 3), n = c(0, 0))) %>%
ungroup() %>%
ggplot() +
geom_bar(aes(x = max_breeding_score, y = n, fill = sex),
alpha = 0.5, stat = "identity", position = "dodge", width = 0.5) +
luke_theme +
theme(legend.position = c(0.15, 0.9)) +
xlab("Maximum individual rufous band score during breeding season") +
ylab("Frequency") +
scale_colour_brewer(palette = "Dark2", direction = -1,
name = "Sex",
labels = c("Female (N = 49)", "Male (N = 35)")) +
scale_fill_brewer(palette = "Dark2", direction = -1,
name = "Sex",
labels = c("Female (N = 49)", "Male (N = 35)"))Males tend to have a more full and immaculate rufous bands than females
effects_mod_max_score <-
as.data.frame(allEffects(mod_max_score)$sex) %>%
mutate(x_mean = ifelse(sex == "F", as.numeric(factor(sex)) + 0.2,
ifelse(sex == "M", as.numeric(factor(sex)) - 0.2, NA)),
x_data = ifelse(sex == "F", as.numeric(factor(sex)),
ifelse(sex == "M", as.numeric(factor(sex)), NA)))
ggplot() +
geom_errorbar(data = effects_mod_max_score,
aes(x = x_mean, ymin = lower, ymax = upper), width = 0.05) +
geom_point(data = effects_mod_max_score,
aes(x = x_mean, y = fit, fill = sex), shape = 21, size = 5) +
geom_jitter(data = ind_breeding_scores,
aes(x = as.numeric(factor(sex)),
y = max_breeding_score, fill = sex),
width = 0.05, height = 0.1, shape = 21, size = 3, alpha = 0.75) +
theme(legend.position = "none",
text = element_text(family = "lato"),
axis.ticks = element_blank(),
axis.text = element_text(size = 10, colour = "grey40"),
axis.title.y = element_text(size = 11, colour = "grey40"),
axis.text.y = element_text(size = 10, colour = "grey40"),
plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
panel.spacing = unit(1, "cm"),
panel.background = element_rect(fill = 'transparent', color = NA),
plot.background = element_rect(fill = 'transparent', color = NA),
strip.background = element_rect(fill = 'grey90', color = NA),
strip.text = element_text(size = 13, colour = "grey40"),
panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.title.x = element_blank()) +
scale_x_continuous(limits = c(0.8, 2.2),
breaks = c(1.1, 1.9),
labels = c("female", "male")) +
ylab("maximum breeding plumage score observed") +
scale_fill_brewer(palette = "Set2") +
scale_y_continuous(limits = c(3.8, 7.2))no effect of age on the maximum extent of an individual’s rufous band score while controlling for underlying sex-specific variation
# draw gt table
mod_max_score_table_status_sex <-
allCoefs_mod %>%
dplyr::select(effect, comp_name, estimate, coefString) %>%
gt(rowname_col = "row",
groupname_col = "effect") %>%
cols_label(comp_name = html("<i>Banded Dotterel rufous band score</i>"),
estimate = "Mean estimate",
coefString = "95% confidence interval") %>%
fmt_number(columns = c(estimate),
rows = 1:10,
decimals = 2,
use_seps = FALSE) %>%
fmt_number(columns = c(estimate),
rows = 11:13,
decimals = 0,
use_seps = FALSE) %>%
sub_missing(columns = 1:4,
missing_text = "") %>%
cols_align(align = "left",
columns = c(comp_name)) %>%
tab_options(row_group.font.weight = "bold",
row_group.background.color = brewer.pal(9,"Greys")[3],
table.font.size = 12,
data_row.padding = 3,
row_group.padding = 4,
summary_row.padding = 2,
column_labels.font.size = 14,
row_group.font.size = 12,
table.width = pct(80))no effect of migratory status on an individual’s rufous band score while controlling for underlying sex-specific variation
# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in sex while accounting for sex differences in maximum rufous band score.
ind_prop_molt_scores <-
dat %>%
mutate(score = as.numeric(score)) %>%
left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>%
filter(!is.na(max_breeding_score)) %>%
mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
select(rings_comb, season, date_J, sex, score, max_breeding_score, prop_molt_score) %>%
arrange(rings_comb, season, date_J) %>%
distinct() %>%
filter(!is.na(prop_molt_score)) %>%
group_by(sex) %>%
mutate(std_max_breeding_score = (max_breeding_score - mean(max_breeding_score, na.rm = TRUE)) / sd(max_breeding_score, na.rm = TRUE)) %>%
ungroup()
# Assess sample sizes of each sex
ind_prop_molt_scores %>%
group_by(sex) %>%
summarise(n_distinct(rings_comb))# A tibble: 2 × 2
sex `n_distinct(rings_comb)`
<chr> <int>
1 F 49
2 M 35
# mixed effects binomial model comparing sex and date effect on the changes in moult scores while controlling for between sex differences in the maximum breeding score
mod1 <-
lme4::glmer(prop_molt_score ~
date_J * sex + std_max_breeding_score +
(1 | rings_comb),
data = ind_prop_molt_scores,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))
# strong date effect, but no sex effect
tbl_regression(mod1, intercept = TRUE,
label = list(date_J ~ "Date", sex ~ "Sex"))| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | 20 | 15, 24 | <0.001 |
| Date | -0.09 | -0.11, -0.07 | <0.001 |
| Sex | |||
| F | — | — | |
| M | 1.1 | -5.7, 8.0 | 0.8 |
| std_max_breeding_score | -0.29 | -0.65, 0.07 | 0.11 |
| Date * Sex | |||
| Date * M | -0.01 | -0.04, 0.03 | 0.7 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
# extract predicted trends
pr <- ggeffects::predict_response(mod1, c("date_J [30:293]", "sex"))
# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date
dates_for_plot <-
data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
date_J = c(1:365))
# join the back-transformed dates to model fits
mod1_fits <-
as.data.frame(pr) %>%
rename(date_J = x,
sex = group) %>%
left_join(., dates_for_plot, by = "date_J")
ind_prop_molt_scores <-
ind_prop_molt_scores %>%
left_join(., dates_for_plot, by = "date_J")
# plot the model
ggplot() +
geom_line(data = mod1_fits,
aes(x = date, y = predicted, color = sex)) +
geom_ribbon(data = mod1_fits,
aes(x = date, ymax = conf.high, ymin = conf.low, fill = sex),
lwd = 1, alpha = 0.25) +
geom_jitter(data = ind_prop_molt_scores,
aes(x = date, y = prop_molt_score, color = sex),
alpha = 0.25, shape = 20, width = 0, height = 0.025) +
theme(legend.position = c(0.3, 0.2),
legend.justification = c(1, 0),
text = element_text(family = "lato"),
axis.ticks = element_blank(),
axis.text.x = element_text(size = 10,
angle = 45,
hjust = 1,
vjust = 1),
axis.text.y = element_text(size = 10, colour = "grey40"),
axis.title.y = element_text(size = 11, colour = "grey40"),
plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
panel.spacing = unit(1, "cm"),
panel.background = element_rect(fill = 'transparent', color = NA),
plot.background = element_rect(fill = 'transparent', color = NA),
panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"),
axis.title.x = element_blank()) +
ylab("proportion intact of individual-based maximum rufous band") +
scale_x_date(date_labels = "%B",
expand = c(0.01, 0.01),
date_breaks = "1 month") +
scale_colour_brewer(palette = "Dark2", direction = -1,
name = "Sex",
labels = c("Female (N = 49)", "Male (N = 35)")) +
scale_fill_brewer(palette = "Dark2", direction = -1,
name = "Sex",
labels = c("Female (N = 49)", "Male (N = 35)"))# wrangle the dates of plumage retention and moult for subsequent analysis
# Step 1: Calculate the max_score for each individual and season
max_scores <- ind_prop_molt_scores %>%
group_by(rings_comb, season) %>%
summarise(max_score = max(prop_molt_score, na.rm = TRUE), .groups = "drop")
# Step 2: Merge the max_scores back into the original data
mid_point_dates <-
ind_prop_molt_scores %>%
left_join(max_scores, by = c("rings_comb", "season")) %>% # Add the max_score
group_by(rings_comb, season) %>%
arrange(date_J) %>% # Ensure data is sorted by date_J
# Step 3: Find the latest date_J for the maximum prop_molt_score
filter(prop_molt_score == max_score) %>%
slice_max(order_by = date_J, n = 1, with_ties = FALSE) %>%
rename(latest_max_date_J = date_J) %>%
# Step 4: Find the earliest date_J where prop_molt_score < max_score
right_join(
ind_prop_molt_scores %>%
left_join(max_scores, by = c("rings_comb", "season")) %>% # Add the max_score here too
group_by(rings_comb, season) %>%
arrange(date_J) %>%
filter(prop_molt_score < max_score) %>%
slice_min(order_by = date_J, n = 1, with_ties = FALSE) %>%
rename(earliest_decline_date_J = date_J),
by = c("rings_comb", "season")
) %>%
# Step 5: Compute the mid-point date_J
mutate(mid_point_date_J = round((latest_max_date_J + earliest_decline_date_J) / 2)) %>%
# Select relevant columns
select(rings_comb, season, latest_max_date_J, earliest_decline_date_J, mid_point_date_J) %>%
distinct()
dates_for_plot <-
data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
date_J = c(1:365))
moult_commence_ring_season <-
left_join(mid_point_dates,
dates_for_plot %>% rename(mid_point_date_J = date_J),
by = "mid_point_date_J")# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in migratory status.
ind_prop_molt_scores <-
dat %>%
mutate(score = as.numeric(score)) %>%
left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>%
filter(!is.na(max_breeding_score)) %>%
mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
select(rings_comb, season, date_J, sex, migratory_status, score, max_breeding_score, prop_molt_score) %>%
arrange(rings_comb, season, date_J) %>%
distinct() %>%
filter(!is.na(prop_molt_score)) %>%
filter(migratory_status %in% c("R", "M"))
# Assess sample sizes of each migratory status
ind_prop_molt_scores %>%
group_by(migratory_status) %>%
summarise(n_distinct(rings_comb))# A tibble: 2 × 2
migratory_status `n_distinct(rings_comb)`
<chr> <int>
1 M 5
2 R 16
# mixed effects binomial model comparing migratory status and date effect on the changes in moult scores
mod1_status <-
lme4::glmer(prop_molt_score ~
date_J * migratory_status +
(1 | rings_comb),
data = ind_prop_molt_scores,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))
# strong date effect, but no status effect
tbl_regression(mod1_status, intercept = TRUE,
label = list(date_J ~ "Date", migratory_status ~ "Migratory status"))| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | 25 | 10, 40 | 0.001 |
| Date | -0.11 | -0.19, -0.04 | 0.002 |
| Migratory status | |||
| M | — | — | |
| R | 12 | -17, 41 | 0.4 |
| Date * Migratory status | |||
| Date * R | -0.06 | -0.20, 0.08 | 0.4 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
# extract predicted trends
pr <- ggeffects::predict_response(mod1_status, c("date_J [30:293]", "migratory_status"))
# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date
dates_for_plot <-
data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
date_J = c(1:365))
# join the back-transformed dates to model fits
mod1_fits <-
as.data.frame(pr) %>%
rename(date_J = x,
migratory_status = group) %>%
left_join(., dates_for_plot, by = "date_J")
ind_prop_molt_scores <-
ind_prop_molt_scores %>%
left_join(., dates_for_plot, by = "date_J")
# plot the model
ggplot() +
geom_line(data = mod1_fits,
aes(x = date, y = predicted, color = migratory_status)) +
geom_ribbon(data = mod1_fits,
aes(x = date, ymax = conf.high, ymin = conf.low, fill = migratory_status),
lwd = 1, alpha = 0.25) +
geom_jitter(data = ind_prop_molt_scores,
aes(x = date, y = prop_molt_score, color = migratory_status),
alpha = 0.25, shape = 20, width = 0, height = 0.025) +
theme(legend.position = c(0.3, 0.2),
legend.justification = c(1, 0),
text = element_text(family = "lato"),
axis.ticks = element_blank(),
axis.text.x = element_text(size = 10,
angle = 45,
hjust = 1,
vjust = 1),
axis.text.y = element_text(size = 10, colour = "grey40"),
axis.title.y = element_text(size = 11, colour = "grey40"),
plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
panel.spacing = unit(1, "cm"),
panel.background = element_rect(fill = 'transparent', color = NA),
plot.background = element_rect(fill = 'transparent', color = NA),
panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"),
axis.title.x = element_blank()) +
ylab("proportion intact of individual-based maximum rufous band") +
scale_x_date(date_labels = "%B",
expand = c(0.01, 0.01),
date_breaks = "1 month") +
scale_colour_brewer(palette = "Set1", direction = -1,
name = "Migratory status",
labels = c("Resident (N = 16)", "Migratory (N = 5)")) +
scale_fill_brewer(palette = "Set1", direction = -1,
name = "Migratory status",
labels = c("Resident (N = 16)", "Migratory (N = 5)"))# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation between age groups.
ind_prop_molt_scores <-
dat %>%
mutate(score = as.numeric(score)) %>%
left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>%
filter(!is.na(max_breeding_score)) %>%
mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
select(rings_comb, season, date_J, sex, age, score, max_breeding_score, prop_molt_score) %>%
arrange(rings_comb, season, date_J) %>%
distinct() %>%
filter(!is.na(prop_molt_score)) %>%
filter(!is.na(age))
# Assess sample sizes of each sex
ind_prop_molt_scores %>%
group_by(age) %>%
summarise(n_distinct(rings_comb))# A tibble: 6 × 2
age `n_distinct(rings_comb)`
<dbl> <int>
1 1 6
2 2 5
3 3 12
4 4 13
5 5 4
6 6 2
# mixed effects binomial model comparing sex and date effect on the changes in moult scores
mod1_age <-
lme4::glmer(prop_molt_score ~
date_J + age +
# date_J * sex +
(1 | rings_comb),# + (1 | season),
data = ind_prop_molt_scores,
family = binomial)#,
# control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))
plot(allEffects(mod1_age))# strong date effect, but no status effect
tbl_regression(mod1_age, intercept = TRUE,
label = list(date_J ~ "Date", age ~ "Age"))| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | 53 | 18, 88 | 0.003 |
| Date | -0.28 | -0.47, -0.10 | 0.003 |
| Age | 2.2 | 0.61, 3.7 | 0.006 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
# extract predicted trends
pr <- ggeffects::predict_response(mod1_age, c("date_J [30:293]", "age"))
# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date
dates_for_plot <-
data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
date_J = c(1:365))
# join the back-transformed dates to model fits
mod1_fits <-
as.data.frame(pr) %>%
rename(date_J = x,
age = group) %>%
left_join(., dates_for_plot, by = "date_J")
ind_prop_molt_scores <-
ind_prop_molt_scores %>%
left_join(., dates_for_plot, by = "date_J")
# plot the model
ggplot() +
geom_line(data = mod1_fits,
aes(x = date, y = predicted, color = age)) +
geom_ribbon(data = mod1_fits,
aes(x = date, ymax = conf.high, ymin = conf.low, fill = age),
lwd = 1, alpha = 0.25) +
geom_jitter(data = ind_prop_molt_scores,
aes(x = date, y = prop_molt_score, color = as.factor(age)),
alpha = 0.25, shape = 20, width = 0, height = 0.025) +
theme(legend.position = c(0.3, 0.2),
legend.justification = c(1, 0),
text = element_text(family = "lato"),
axis.ticks = element_blank(),
axis.text.x = element_text(size = 10,
angle = 45,
hjust = 1,
vjust = 1),
axis.text.y = element_text(size = 10, colour = "grey40"),
axis.title.y = element_text(size = 11, colour = "grey40"),
plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
panel.spacing = unit(1, "cm"),
panel.background = element_rect(fill = 'transparent', color = NA),
plot.background = element_rect(fill = 'transparent', color = NA),
panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"),
axis.title.x = element_blank()) +
ylab("proportion intact of individual-based maximum rufous band") +
scale_x_date(date_labels = "%B",
expand = c(0.01, 0.01),
date_breaks = "1 month",
limits = c(as.Date("2021-10-01"), as.Date("2022-05-01"))) +
scale_colour_brewer(palette = "Spectral", direction = -1,
name = "Age",
labels = c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)",
"4 (N = 13", "5 (N = 4)", "6 (N = 2)")) +
scale_fill_brewer(palette = "Spectral", direction = -1,
name = "Age",
labels = c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)",
"4 (N = 13", "5 (N = 4)", "6 (N = 2)"))# explore datasets provided by Ted
## 2021 season
# import and consolidate key columns
breeding_data_2021 <-
read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_2021_Feb2022_upd_Nov22_Locations.xlsx",
sheet = "Kaikoura_NestVisit2021_1",
col_types = "text") %>%
mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>%
mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>%
mutate(date = as.Date(visit_date_nz)) %>%
select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>%
rename(nest_fate = `Nest Status`,
eggs = `Egg count`,
hatched = `Number Hatched`,
fledged = `Number Fledged`,
chicks = `Number of chicks seen`)
# Define a function to extract the desired patterns and ensure enough columns
extract_patterns <- function(text) {
# Replace vertical bar '|' with an empty string
text <- gsub("\\|", "", text)
# Extract all 4 or 5-character long texts that start with R or r
matches_r <- str_extract_all(text, "\\b[Rr]\\w{3,4}\\b")[[1]]
# Extract UB, UN, UBF, UBM as standalone patterns
matches_ub_un <- str_extract_all(text, "\\b[Uu][BbNnFfMm]\\b")[[1]]
# Combine matches
matches <- c(matches_r, matches_ub_un)
# Ensure all matches are capitalized
matches <- toupper(matches)
# Ensure there are exactly two columns, filling with NA if necessary
result <- c(matches, rep(NA, 2 - length(matches)))
# Return the result as a named vector
return(setNames(result, c("parent1", "parent2")))
}
# Apply the function to the text column and store the results in new columns
extracted_parents <- t(sapply(breeding_data_2021$NestID, extract_patterns))
breeding_data_2021 <-
bind_cols(breeding_data_2021, extracted_parents) %>%
select(-parent2) %>%
rename(parent = parent1) %>%
bind_rows(bind_cols(breeding_data_2021, extracted_parents) %>%
select(-parent1) %>%
rename(parent = parent2)) %>%
filter(!is.na(parent)) %>%
filter(parent %in% (dat %>% filter(season == 2021) %>% pull(rings_comb) %>% unique())) %>%
filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>%
group_by(parent) %>%
arrange(desc(date)) %>%
slice(1) %>%
ungroup() %>%
arrange(desc(date)) %>%
mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>%
rename(last_date_breeding = date) %>%
select(parent, last_date_breeding, nest_id, OBJECTID)
# import and consolidate key columns
breeding_data_2021_ <-
read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_2021_Feb2022_upd_Nov22_Locations.xlsx",
sheet = "Kaikoura_DotterelNest2021_0",
col_types = "text") %>%
mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"),
fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"),
found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>%
mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"),
fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"),
found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>%
mutate(hatch_date = as.Date(hatch_date_nz),
fail_date = as.Date(fail_date_nz),
found_date = as.Date(found_date_nz)) %>%
mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>%
select(OBJECTID, nest_id, found_date, hatch_date, fail_date)
breeding_2021_final <- left_join(breeding_data_2021, breeding_data_2021_, by = "nest_id") %>%
mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) |
((!is.na(hatch_date) & found_date < hatch_date)) |
((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) |
((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>%
group_by(parent) %>%
mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>%
mutate(days_diff = last_date_breeding_ - last_date_breeding) %>%
arrange(desc(days_diff)) %>%
mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding))) %>%
arrange(desc(last_date_breeding_final)) %>%
select(parent, last_date_breeding_final) %>%
# specify the season as the first calender year
mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final)))
###########
## 2022 season
# import and consolidate key columns
breeding_data_2022 <-
read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Jan23_upd.xlsx",
sheet = "Kaikoura_NestVisit_1",
col_types = "text") %>%
mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>%
mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>%
mutate(date = as.Date(visit_date_nz)) %>%
select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>%
rename(nest_fate = `Nest Status`,
eggs = `Egg count`,
hatched = `Number Hatched`,
fledged = `Number Fledged`,
chicks = `Number of chicks seen`)
# Apply the function to the text column and store the results in new columns
extracted_parents2 <- t(sapply(breeding_data_2022$NestID, extract_patterns))
breeding_data_2022 <-
bind_cols(breeding_data_2022, extracted_parents2) %>%
select(-parent2) %>%
rename(parent = parent1) %>%
bind_rows(bind_cols(breeding_data_2022, extracted_parents2) %>%
select(-parent1) %>%
rename(parent = parent2)) %>%
filter(!is.na(parent)) %>%
filter(parent %in% (dat %>% filter(season == 2022) %>% pull(rings_comb) %>% unique())) %>%
filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>%
group_by(parent) %>%
arrange(desc(date)) %>%
slice(1) %>%
ungroup() %>%
arrange(desc(date)) %>%
mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>%
rename(last_date_breeding = date) %>%
select(parent, last_date_breeding, nest_id, OBJECTID)
# import and consolidate key columns
breeding_data_2022_ <-
read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Jan23_upd.xlsx",
sheet = "Kaikoura_DotterelNest_0",
col_types = "text") %>%
mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"),
fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"),
found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>%
mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"),
fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"),
found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>%
mutate(hatch_date = as.Date(hatch_date_nz),
fail_date = as.Date(fail_date_nz),
found_date = as.Date(found_date_nz)) %>%
mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>%
select(OBJECTID, nest_id, found_date, hatch_date, fail_date)
breeding_2022_final <- left_join(breeding_data_2022, breeding_data_2022_, by = "nest_id") %>%
mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) |
((!is.na(hatch_date) & found_date < hatch_date)) |
((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) |
((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>%
group_by(parent) %>%
mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>%
mutate(days_diff = last_date_breeding_ - last_date_breeding) %>%
arrange(desc(days_diff)) %>%
mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding))) %>%
arrange(desc(last_date_breeding_final)) %>%
select(parent, last_date_breeding_final) %>%
# specify the season as the first calender year
mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final)))
###########
## 2023 season
# import and consolidate key columns
breeding_data_2023 <-
read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Feb2024_upd.xlsx",
sheet = "Kaikoura_NestVisit_1",
col_types = "text") %>%
mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>%
mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>%
mutate(date = as.Date(visit_date_nz)) %>%
select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>%
rename(nest_fate = `Nest Status`,
eggs = `Egg count`,
hatched = `Number Hatched`,
fledged = `Number Fledged`,
chicks = `Number of chicks seen`)
# Apply the function to the text column and store the results in new columns
extracted_parents3 <- t(sapply(breeding_data_2023$NestID, extract_patterns))
breeding_data_2023 <-
bind_cols(breeding_data_2023, extracted_parents3) %>%
select(-parent2) %>%
rename(parent = parent1) %>%
bind_rows(bind_cols(breeding_data_2023, extracted_parents3) %>%
select(-parent1) %>%
rename(parent = parent2)) %>%
filter(!is.na(parent)) %>%
filter(parent %in% (dat %>% filter(season == 2023) %>% pull(rings_comb) %>% unique())) %>%
filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>%
group_by(parent) %>%
arrange(desc(date)) %>%
slice(1) %>%
ungroup() %>%
arrange(desc(date)) %>%
mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>%
rename(last_date_breeding = date) %>%
select(parent, last_date_breeding, nest_id, OBJECTID)
# import and consolidate key columns
breeding_data_2023_ <-
read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Feb2024_upd.xlsx",
sheet = "Kaikoura_DotterelNest_0",
col_types = "text") %>%
mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"),
fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"),
found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>%
mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"),
fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"),
found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>%
mutate(hatch_date = as.Date(hatch_date_nz),
fail_date = as.Date(fail_date_nz),
found_date = as.Date(found_date_nz)) %>%
mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>%
select(OBJECTID, nest_id, found_date, hatch_date, fail_date)
breeding_2023_final <- left_join(breeding_data_2023, breeding_data_2023_, by = "nest_id") %>%
mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) |
((!is.na(hatch_date) & found_date < hatch_date)) |
((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) |
((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>%
group_by(parent) %>%
mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>%
mutate(days_diff = last_date_breeding_ - last_date_breeding) %>%
arrange(desc(days_diff)) %>%
mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding))) %>%
arrange(desc(last_date_breeding_final)) %>%
select(parent, last_date_breeding_final) %>%
# specify the season as the first calender year
mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final)))
# combine breeding data from the 3 seasons into one table
breeding_data_all <- bind_rows(breeding_2021_final,breeding_2022_final,breeding_2023_final) %>%
rename(rings_comb = parent) %>%
# change to Julian date shifted for the Southern Hemisphere (1 = July 1)
mutate(breeding_date_J = as.numeric(format(last_date_breeding_final + 181, "%j"))) %>%
ungroup() %>%
# standardize the last breeding date by year
mutate(breeding_date_J_s = scale_by(breeding_date_J ~ season, ., scale = 0)) %>%
# make the scaled date variable numeric class
mutate(breeding_date_J_snum = as.numeric(breeding_date_J_s))ind_prop_molt_scores_breeding_date <-
dat %>%
mutate(score = as.numeric(score)) %>%
left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>%
filter(!is.na(max_breeding_score)) %>%
mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
select(rings_comb, season, date_J, sex, migratory_status, age, score, max_breeding_score, prop_molt_score) %>%
arrange(rings_comb, season, date_J) %>%
distinct() %>%
filter(!is.na(prop_molt_score)) %>%
full_join(breeding_data_all, by = c("rings_comb", "season")) %>%
mutate(date_diff = abs(breeding_date_J - date_J)) %>%
arrange(rings_comb, season, date_diff) %>%
group_by(rings_comb, season) %>% slice(1) %>% filter(date_diff < 30) %>%
dplyr::select(rings_comb, season, sex, breeding_date_J, breeding_date_J_snum, prop_molt_score, max_breeding_score, age) %>%
distinct() %>%
group_by(sex) %>%
mutate(std_max_breeding_score = (max_breeding_score - mean(max_breeding_score, na.rm = TRUE)) / sd(max_breeding_score, na.rm = TRUE)) %>%
ungroup() %>%
mutate(breeding_date_J_s = scale(breeding_date_J)) %>%
left_join(., moult_commence_ring_season, by = c("rings_comb", "season"))do individuals with a fuller breeding plumage (for their sex) breed later than the rest of the population? No effect
mod_max_score_breeding_date <-
lme4::lmer(breeding_date_J_snum ~
std_max_breeding_score +
(1 | rings_comb),
data = ind_prop_molt_scores_breeding_date)
tbl_regression(mod_max_score_breeding_date, intercept = TRUE)| Characteristic | Beta | 95% CI1 |
|---|---|---|
| (Intercept) | 4.5 | -0.64, 9.7 |
| std_max_breeding_score | 2.5 | -2.7, 7.7 |
| 1 CI = Confidence Interval | ||
plot(allEffects(mod_max_score_breeding_date), type = "response")do individuals with a late breeding attempts tend to have a fuller breeding plumage? No effect
mod_breeding_date_max_score <-
lme4::lmer(std_max_breeding_score ~
breeding_date_J_snum +
(1 | rings_comb),
data = ind_prop_molt_scores_breeding_date)
summary(mod_max_score_breeding_date)Linear mixed model fit by REML ['lmerMod']
Formula: breeding_date_J_snum ~ std_max_breeding_score + (1 | rings_comb)
Data: ind_prop_molt_scores_breeding_date
REML criterion at convergence: 1075.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.7013 -0.6280 0.1192 0.6279 1.9243
Random effects:
Groups Name Variance Std.Dev.
rings_comb (Intercept) 0.0 0.0
Residual 795.4 28.2
Number of obs: 114, groups: rings_comb, 63
Fixed effects:
Estimate Std. Error t value
(Intercept) 4.539 2.641 1.718
std_max_breeding_score 2.475 2.665 0.929
Correlation of Fixed Effects:
(Intr)
std_mx_brd_ 0.000
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
plot(allEffects(mod_breeding_date_max_score), type = "response")do individuals with late breeding attempts have a delayed moult? Yes: individuals that breed later than the rest of the population in a given year tend to commence their complete/post-breeding moult later
mod_breeding_date_prop_molt_mid <-
lme4::lmer(mid_point_date_J ~
breeding_date_J_snum +
(1 | rings_comb),
data = ind_prop_molt_scores_breeding_date)
tbl_regression(mod_breeding_date_prop_molt_mid, intercept = TRUE)| Characteristic | Beta | 95% CI1 |
|---|---|---|
| (Intercept) | 153 | 145, 160 |
| breeding_date_J_snum | 0.31 | 0.03, 0.58 |
| 1 CI = Confidence Interval | ||
plot(allEffects(mod_breeding_date_prop_molt_mid), type = "response")do individuals with late breeding attempts retain their full breeding plumage longer? Yes: individuals that breed later than the rest of the population in a given year tend to retain their full breeding plumage longer
mod_breeding_date_prop_molt_latest <-
lm(latest_max_date_J ~ breeding_date_J_snum,
data = ind_prop_molt_scores_breeding_date)
tbl_regression(mod_breeding_date_prop_molt_latest, intercept = TRUE)| Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | 148 | 141, 155 | <0.001 |
| breeding_date_J_snum | 0.59 | 0.34, 0.85 | <0.001 |
| 1 CI = Confidence Interval | |||
plot(allEffects(mod_breeding_date_prop_molt_latest), type = "response")# wrangle data to get ages and dates of last breeding for each season of all
# birds available
age_breeding_date <-
dat %>%
select(rings_comb, season, age, sex) %>%
full_join(breeding_data_all, by = c("rings_comb", "season")) %>%
distinct()do older individuals conclude their breeding season later? No: there is strong evidence that younger individuals tend to have the latest breeding attempts in the population for a given year.
mod_age_breeding_date <-
lme4::lmer(breeding_date_J_snum ~ age +
(1 | rings_comb),
data = age_breeding_date)
summary(mod_age_breeding_date)Linear mixed model fit by REML ['lmerMod']
Formula: breeding_date_J_snum ~ age + (1 | rings_comb)
Data: age_breeding_date
REML criterion at convergence: 312.2
Scaled residuals:
Min 1Q Median 3Q Max
-2.3355 -0.5314 0.1108 0.5764 2.2836
Random effects:
Groups Name Variance Std.Dev.
rings_comb (Intercept) 0 0.0
Residual 1069 32.7
Number of obs: 33, groups: rings_comb, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 25.976 14.414 1.802
age -7.608 3.468 -2.194
Correlation of Fixed Effects:
(Intr)
age -0.919
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
plot(allEffects(mod_age_breeding_date), type = "response")# strong date effect, but no status effect
tbl_regression(mod_age_breeding_date, intercept = TRUE,
label = list(age ~ "Age"))| Characteristic | Beta | 95% CI1 |
|---|---|---|
| (Intercept) | 26 | -2.3, 54 |
| Age | -7.6 | -14, -0.81 |
| 1 CI = Confidence Interval | ||
# extract predicted trends
# I tried the whole range for the breeding timing (breeding_date_J [72:221]), but then decided to take some values with 20 days interval to visualize it better
pr <- ggeffects::predict_response(mod_age_breeding_date, c("age"))
dates_for_plot <-
data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
breeding_date_ = c(1:365))
CL_dates_for_plot <-
data.frame(conf.low_date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
conf.low_ = c(1:365))
CH_dates_for_plot <-
data.frame(conf.high_date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
conf.high_ = c(1:365))
# join the back-transformed dates to model fits
mod1_fits <-
as.data.frame(pr) %>%
rename(age = x,
# sex = group,
breeding_date = predicted) %>%
mutate(breeding_date_ = round(breeding_date),
conf.low_ = round(conf.low),
conf.high_ = round(conf.high)) %>%
left_join(., dates_for_plot, by = "breeding_date_") %>%
left_join(., CL_dates_for_plot, by = "conf.low_") %>%
left_join(., CH_dates_for_plot, by = "conf.high_")
# plot the model
ggplot() +
geom_line(data = mod1_fits,
aes(x = age, y = date)) +
geom_ribbon(data = mod1_fits,
aes(x = age, ymax = conf.high_date, ymin = conf.low_date),
lwd = 1, alpha = 0.25) +
luke_theme +
# theme(legend.position = c(0.3, 0.2),
# legend.justification = c(1, 0),
# strip.background = element_blank(),
# axis.title.x = element_blank(),
# axis.text.x = element_text(size = 10,
# angle = 45,
# hjust = 1,
# vjust = 1)) +
ylab("date of last breeding")add the breeding data to the analysis: for each individual include the date of last breeding (i.e., last date seen with nest or brood) as a predictor of molt timing. Also do an analysis relating the timing of last breeding to the age of the individual. Bashar will send Luke the breeding data he received from Ted.
drop pre-October data from the molt timing analysis.
next steps: Bashar to eventually write up this study as a report (which could eventually be used as his MSc thesis if he wants). Priority will be on him writing up a proposal for his thesis (which has a concrete deadline and is graded). Bashar will present the study at our weekly Monday meeting in the first few months of 2025.